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Based on the quasiparticle model of the quark-gluon plasma (QGP), a color quantum path-integral 
Monte-Carlo (PIMC) method for calculation of thermodynamic properties and - closely related to 
the latter - a Wigner dynamics method for calculation of transport properties of the QGP are 
formulated. The QGP partition function is presented in the form of a color path integral with 
a new relativistic measure instead of the Gaussian one traditionally used in the Feynman- Wiener 

■ path integral. A procedure of sampling color variables according to the SU(3) group Haar measure 
is developed for integration over the color variable. It is shown that the PIMC method is able to 

, reproduce the lattice QCD equation of state at zero baryon chemical potential at realistic model 

■ parameters (i.e. quasiparticle masses and coupling constant) and also yields valuable insight into the 
internal structure of the QGP. Our results indicate that the QGP reveals quantum liquid-like (rather 

^ , than gas-like) properties up to the highest considered temperature of 525 MeV. The pair distribution 

■ functions clearly reflect the existence of gluon-gluon bound states, i.e. glueballs, at temperatures 
just above the phase transition, while meson-like qq bound states are not found. The calculated self- 

I diffusion coefficient agrees well with some estimates of the heavy-quark diffusion constant available 

from recent lattice data and also with an analysis of heavy-quark quenching in experiments on 
ultrarelativistic heavy ion collisions, however, appreciably exceeds other estimates. The lattice and 
heavy-quark-quenching results on the heavy-quark diffusion are still rather diverse. The obtained 
results for the shear viscosity are in the range of those deduced from an analysis of the experimental 
elliptic flow in ultrarelativistic heavy ions collisions, i.e. in terms the viscosity-to-entropy ratio, 
O ■ 1/4-K < ri/S < 2.5/4-K, in the temperature range from 170 to 440 MeV. 
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Determining the properties of the quark-gluon plasma (QGP) is nowadays one of the most important goals in high- 
energy nuclear physics. In recent years, experiments at the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven 
. National Laboratory [l| and the Large Hadron Collider (LHC) at CERN have provided a wealth of data from which 

■ one can obtain information on a number of features of the QGP. The most striking result, obtained from analysis of 
^ . these experimental data [U, 01 1 is that the deconfined quark-gluon matter behaves as an almost perfect fluid rather 

• ■ than a perfect gas, as it could be expected from the asymptotic freedom. 

! There are various approaches for a theoretical study of the QGP each of which has its advantages and disadvantages. 
^ ■ The most fundamental way to compute the properties of strongly interacting matter is provided by the lattice QCD 
" " ' In particular, the lattice QCD computations 0] indicated that the transition into the QGP phase at small 

chemical potentials has a nature of a rapid crossover, which, strictly speaking, is not a phase transition. However, 
following the custom of the heavy-ion society, we call this transformation as a phase transition of the crossover 
type. Interpretation of these very complicated computations requires application of various QCD motivated, albeit 
schematic, models simulating various aspects of the full theory. Moreover, such models are needed in cases when the 
lattice QCD fails, e.g. at large quark chemical potentials and out of equilibrium. While certain progress has been 
achieved in recent years, transport properties of the QGP are still poorly accessible within lattice QCD - only viscosity 
for pure gluonic matter was computed Q. It is, therefore, crucial to devise reliable and manageable theoretical tools 
for a quantitative description of non-Abelian QGP both in and out of equilibrium. 

A semi-classical approximation, based on a quasiparticle picture has been introduced in Refs. (9l-[l3j. It is motivated 
by the expectation that the main features of non-Abelian plasmas can be understood in simple semi-classical terms 
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without the difficulties inherent to a full quantum field-theoretical analysis. Independently the same ideas were 
implemented in terms of molecular dynamics (MD) [l^ : this approach was further developed in a series of works 
[Hill. The MD allows one to treat soft processes in the QGP which are not accessible by perturbative means. 

Strongly correlated behavior of the QGP is expected to show up in long-ranged spatial correlations of quarks and 
giuons which, in fact, may give rise to liquid-like and, possibly, solid-like structures. This expectation is based on a 
very similar behavior observed in electrodynamic plasmas |15l llTI . Il8| . This similarity was exploited to formulate a 
classical nonrelativistic model of a color Coulomb interacting QGP |15| which was numerically analyzed by classical 
MD simulations. Quantum effects were either neglected or incorporated phcnomenologically via a short-range repulsive 
correction to the pair potential. Such a rough model may, however, become a critical issue at high densities. Similar 
models in electrodynamic plasmas showed poor behavior in the region of strong wave function overlap, in particular 
around the Mott density where bound states break up. For temperatures and densities of the QGP considered in 
Ref. [l^ these effects are very important since the quasiparticle thermal wave length is of the order of the average 
interparticle distance. Therefore, to account for quantum effects, we follow an idea of Kelbg [l^ that allows one 
to rigorously include quantum corrections to the pair potential^. Strictly speaking, this method is applicable only 
to weak coupling. To extend the method to stronger couplings, an "improved Kelbg potential" was derived, which 
contains a single free parameter, being fitted to the exact solution of the quantum-mechanical two-body problem. 
Using the method of the improved Kelbg potential in classical MD simulations one is able to describe thermodynamic 
properties of a partially ionized plasma up to moderate couplings [2l| . However, this approach may fail, if bound 
states of more than two particles are formed in the system. This is a result of break-down of the pair approximation 



for the density matrix, as demonstrated in Refs. [21| . A superior approach, which does not have this limitation, 
consists in the use of the original Kelbg potential in path integral Monte Carlo simulations (PIMC) which effectively 
map the problem onto a high-temperature weakly coupled and weakly degenerate one. This allows one to advance 
the analysis to strong couplings and is, therefore, a relevant choice for the present purpose. The PIMC method 
has been successfully applied to various phases of strongly coupled electrodynamic plasmas [12, [23l. Examples are 
partially- ionized dense hydrogen plasmas, where liquid-like and crystalline behavior was observed |24| - [26{ . as well as 
electron- hole plasmas in semiconductors [13, H^, including excitonic bound states. 

In this paper we extend the previous classical nonrelativistic simulations [l^ in two ways: first, we include quantum 
and spin effects and, second, we take into account the dominant relativistic effects, cf. section |TT1 This is done in the 
frame of quantum Monte Carlo simulations where we rewrite the partition function of this system in the form of color 
path integrals with a new relativistic measure instead of Gaussian one used in Fcynman- Wiener path integrals. For 
the integration of the partition function over color variables we develop a procedure of sampling the color quasiparticle 
variables according to the SU(3) group Haar measure with the quadratic and cubic Casimir conditions. The developed 
approach self-consistently takes into account the Fermi (Bose) statistics of quarks (giuons). The main goal of this 
article is to test the present Color Path-Integral Monte-Carlo (PIMC) against known lattice data Q and to predict 
additional properties of the QGP, which are still inaccesible from lattice QCD. First results of the path integral 
approach for thermodynamic properties of the nonideal QGP have been already briefly reported in Ref. [231 for SU(3) 
group and in Refs. [3(]| - l33j for SU(2) group. In this paper we show that the PIMC method is able to reproduce the 
QCD lattice equation of state at vanishing baryon-charge density and also yields valuable insight into the internal 
structure of the QGP. These results are presented in section HV Al 

Hydrodynamic simulations of relativistic heavy-ion collisions require not only knowledge of thermodynamic prop- 
erties of the QGP but also of the transport properties. Unfortunately the PIMC method itself is not able to directly 
predict transport properties. Therefore, to simulate quantum QGP transport and thermodynamic properties within 
a unified approach we combine the path integral and Wigner (in phase space) formulations of quantum mechanics 
f section lllip . There the kinetic coefficients are calculated by means of Kubo formulas. In this approach the PIMC 
method is used to generate initial conditions (equilibrium quasiparticle configurations) for dynamical trajectories de- 
scribing the time evolution for spatial, momentum and color variables. Correlation functions and kinetic coefficients 
are calculated as averages of Weyl's symbols of dynamic operators along these trajectories. The basic ideas of this 
approach have been published in Ref. [SJI • Using this approach we calculate the self-diffusion coefficient and viscosity 
of the strongly coupled QGP. These results are presented in section IIVBI 



The idea to use a Kelbg-type effective potential also for quark matter was independently proposed by K. Dusling and C. Young 120 
However, their potentials are limited to weakly nonideal systems. 
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II. THERMODYNAMICS OF QGP 

111 this section we summarize the main ideas of our approach to the thermodynamic properties of the strongly 
correlated quark-gluon plasma. This approach is based on a generalization of the Feynman path integral representation 
of quantum mechanics to high energy matter. Before deriving the main equations of our PIMC approach we specify 
the simplifications and model parameters. 

A. Basics of the model 

The basic assumptions of the model are similar to those of Ref. [3l : 

I. Quasiparticles masses (m) are of order or higher than the mean kinetic energy per particle. This assumption is 

based on the analysis of QCD lattice data [35l - [37j . For instance, at zero net-baryon density it amounts to m > T, 
where T is a temperature. 

II. In view of the first assumption, interparticle interaction is dominated by a color-electric Coulomb potential. 

Magnetic effects are neglected as subleading ones. 

III. Relying on the fact that the color representations are large, the color operators are substituted by their average 

values, i.e. by Wong's classical color vectors (8D in SU(3)) with the quadratic and cubic Casimir conditions 

IV. We consider the 3-flavor quark model. Since the masses of the 'wp', 'down' and 'strange' quarks extracted from 

lattice data are still indistinguishable, we assume these masses to be equal. As for the gluon quasiparticles, we 
allow their mass to be different (heavier) from that of quarks. 

The quality of these approximations and their limitations were discussed in Ref. [T5| . Thus, this model requires the 
following quantities as a function of temperature (T) and quark chemical potential (fiq) as an input: 

1. quasiparticle masses, for quarks mq and gluons nig, and 

2. the coupling constant g^, or as = /An. 

Input quantities should be deduced from lattice QCD data or from an appropriate model simulating these data. 

It has been established that hard modes (in terms of hard thermal loop approximation (s^-jHI) behave like quasi- 
particles (4l| . Therefore, masses of these quasiparticles should be deduced from nonperturbative calculations taking 
into account hard field modes, e.g., they can be associated with pole masses deduced from lattice QCD calculations. 
At the same time, the soft quantum fields are characterized by very high occupation numbers per mode. Therefore, 
to leading order, they can be well approximated by soft classical fields. This is precisely the picture we are going 
to utilize: massive quantum quasiparticles (hard modes) interacting via classical color fields. Applicability of such 
approach was discussed in Refs. [l^, [13 i^i detail. Our approach differs from that of Ref. EBl by a quantum 
treatment to quasiparticles instead of the classical one, and additionally by a relativistic desciption of the kinetic 
energy instead of the nonrelativistic approximation of Ref. [loj . 

B. Color Path Integrals 

We consider a multi-component QGP consisting of N color quasiparticles: Ng gluons, Nq quarks and Nq antiquarks. 
The Hamiltonian of this system is H = K + U'-^ with the kinetic and color Coulomb interaction parts 

Here i and j summations run over quark and gluon quasiparticles, i,j — l,...,iV, N = Nq + Nq + Ng, Nq = 
Nu + Nd. + Ns and Nq = Nu + Nd + Ng are total numbers of quarks and antiquarks of all flavours {up, down and 
strange), respectively, 3D vectors are quasiparticle spatial coordinates, the Qi denote the Wong's quasiparticle color 
variable (8D-vector in the group SU{3)), {Qi-Qj) denote scalar product of color vectors. Nonrelativistic approximation 
for potential energy is used, while for kinetic energy we still keep relativistic form as the quasiparticle masses are not 
negligible as compared with temperature. The eigenvalue equation of this Hamiltonian is usually called the spinless 
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Salpeter equation. It may be regarded as a well-defined approximation to the BetheSalpeter formalism [42| for the 
description of bound states within relativistic quantum field theories, obtained when assuming that all bound-state 
constituents interact instantaneously and propagate like free particles [i^. Among others, it yields semirelativistic 
descriptions of hadrons as bound states of quarks [13, EEl • 

In the classic approximation this system is governed by Wong's equations of motion (ssj 



dpi jt) 

dt 
dr,{t) 

dt 
dQfjt) 

dt 



F.W, (2) 

Mt), (3) 

(4) 



where Pi is the momentum of a quasiparticle, — Pi / i/pf + mf is its velocity, = —dU^/dri is the color-electric 
force experienced by the quasiparticle, and 



is the driving force in equation of motion for the color charge, jabc are structure constants of the group SU(3) and 
a, 6, c = 1, . . . , 8 (see Appendix I). Below we consider spatial degrees of freedom quantum-mechanically while the color 
dynamics, still classically. 

Thermodynamic properties in the grand canonical ensemble with given temperature (T), net-quark- number (/i^) 
and strange (/i^) chemical potentials, and fixed volume V are fully described by the grand partition function^ 

= V exp{,,(iV, - N,)in_ exp{,.^A.. - iV.)/T} ^ 

Z(XN},V,fi) = E / dr diiQ p{r,Q.o-ANY^). (7) 

where {A^} = {N^, N^i, Ng, Nu, Nd, Ns, Ng} , p{r,Q,a;{N}; f3) denotes the diagonal matrix elements of the density- 
matrix operator p = exp(— with /3 = 1/T. Here r, a and Q denote the multi-dimensional vectors related 
spatial, spin and color degrees of freedom, respectively, of all quarks, antiquarks and gluons. The a summation and 
spacial {dr = d^ri...d^rN) and color (d/xQ = dpQi...dpQN) integrations run over all individual degrees of freedom of 
the particles, dpQi denotes integration over SU(3) group Haar measure, see Appendix I. Usual choice of the strange 
chemical potential is ps = —pg (nonstrange matter), such that the total factor in front of {Ns~ Ns) is zero. Therefore, 
below we omit ps from the list of variables. 

Since the masses and the coupling constant depend on the temperature and quark chemical potential, special care 
should be taken to preserve thermodynamical consistency of this approach. To achieve this, thermodynamic functions 
such as pressure, P, entropy, S, baryon number, Nb, and internal energy, E, should be calculated through respective 
derivatives of the logarithm of the partition function 

P = d{T\nZ)/dV, (8) 

S = diT\nZ)/dT, (9) 

Nb = {l/3)diTlnZ)/dpq, (10) 

E = -PV + TS + 3pqNB. (11) 

This is a conventional way of maintaining the thermodynamical consistency in approaches of the Ginzburg-Landau 
type as they are applied in high-energy physics, e.g., in the PNJL model. 



^ Here we explicitly write sum over different quark flavors (u,d,s). Below the sum over quark degrees of freedom is understood in the 
same way. 
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The exact density matrix p = e of interacting quantum systems can be constructed using a path integral 
approach [i^ [4^ based on the operator identity e~^^ = e~^^-^ • e~^^^ . . . e~^^^ , where the r.h.s. contains n + 1 
identical factors with = (3/{n+ 1), which ahows us to rewrite'^ the integral in Eq. ([7]) in the form 

^ / dr(0)dAiQp(r("\Q,a;{7V};/3) 

d^iQ [ dr(0)dr(i) ...dr("+i)i?(r(0),r(i\...,r("+i);Q;{iV};/3) (5 (r("+i) -r(°)) . (12) 



Notice that the color charge Q is a classical variable already in the mixed (i.e. coordinate-momentum) representation, 
see Appendix I. Therefore, we do not build a chain of n different Q-variables. The spin variable a is the same in all 
p^'-* except for the p'""*"^^, where it is initially set cr*-""'"^^ and only after permutations performed is put to a. The spin 
gives rise to the spin part of the density matrix (S) with exchange effects accounted for by the permutation operators 
Pq, Pq and Pg acting on spatial A"-'^'^'^ degrees of freedom and the spin projections cr'^"'"'""'^^. The sum runs over all 
permutations with parity Kp^ and Kp_. In Eq. ([T2|) 



P , rW , Q; {TV}; A/?) = (rC-^ , Q je"^^* | r« , q) , (13) 



is an off-diagonal element of the density matrix. Accordingly each quasiparticle is represented by a set of coordi- 
nates {r^P , . . . ,7'-"^} ("beads") and a 8-dimensional color vector Qi in the SU{?>) group. Thus, all "beads" of each 
quasiparticle are characterized by the same spin projection, flavor and color charge. Notice that masses and coupling 
constant in each p^'^ are the same as those for the original quasiparticles, i.e. these are still defined by the actual 
temperature T. 

The main advantage of decomposition (|12p is that it allows us to use perturbation theory to obtain approximation 
for density matrices p''', which is applicable due to smallness of artificially introduced factor l/(n + 1). This means 
that in each p(') the ratio g'^{T, fiq){Q, ■ Qj)/[47r|rf ^ - rj'' \T{n + l)] can be always made much smaller than one, which 
allows us to use perturbation theory with respect to the potential. Each factor in Eq. p2|) should be calculated with 
the accuracy of order of l/(ri, -I- 1)^ with > 1, as in this case the error of the whole product in the limit of large n will 
be equal to zero. In the limit {n + 1) — > oo, p' can be approximated by a product of two-particle density matrices 
pfj [IE Si 113 ■ This approximation can be deduces from operator expansion 

exp (^-A^H^ = exp (^~A/3K^ exp (^-A^[/^) exp (-A/?^ K, /2^... (14) 

As the first approximation with the error proportional to (l/(n + 1))^ , we can write 



exp 



-A^ij) « exp (^-Apfdj exp (^-Apu'^^ , (15) 



It is very important that in this case the error of the whole product in Eq. (jl2p is proportianal to the l/(ri + 1) and is 
equal to zero in the limit n — > oo. The second advantage of the decomposition in Eq. (jl2p is that it reduces quantum 
multi-particle interaction to the pair-wise sum of two-particle interactions described by two-particle classical density 
matrices in each p^'^. Thus, neglecting the commutator terms in Eq. (|14p. we arrive at the following expression for 
the density matrix of Eq. ((T3)) 



(0 _ (0 

P neglecting commutators ^0 *--^P 



A/3C/'^ (rW,Q)] , (16) 



where pg ■* is the corresponding density matrix of noninteracting particles. This approximation works well for potentials 
bounded below. However, the Coulomb potential can go to minus infinity and hence the result (jl6p diverges in this 
limit. 



^ For the sake of notation convenience, we ascribe superscript to the original variables. 
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A more sophisticated treatment is required to avoid this divergence. Ah the calculations along this line can be 
rigorously performed for the two-particle density matrix pl^l (r, r', Q; A/3), where r = {ri,r2}, r' = {r'j^jr'j} and 
Q = {Qi-,Q2}- Expanding the two-particle density matrix up to the second order in l/{n + 1), one arrives to the 
following result [l^ 



p[2](r,r',Q;A/3)s,p[2l(r,r',g;A/3) 



dr / d-'r- 



47r|?|AA22v/7(W 



: exp 



7r|ri2 - r|" 



exp 



7r|?-r'i2p 



pjfl exp hA/3 $i2(ri2, r'i2, Qi,Q2)] 



(17) 



where ri2 = ri — r2, r'12 = r'l — r'2, AA12 = ■\/27rA/3/mi2 is defined in terms of the reduced mass of the pair of 

[2] 

particles: mi2 — mim2 / {mi + TO2), and Pq is the two-particle density matrix of noninteracting particles. In the end 
of Eq. (|17p the result is presented in the form similar to E g. (iT5|) . i.e. in terms of an off-diagonal two-particle effective 
quantum potential $12, which is called a Kelbg potential Eq. ((T7| is the definition of the color Kelbg potential. 
The diagonal part of the color Kelbg potential can be obtained analytically 

9' {Q 



$i2(ri2,ri2,(3i,<92) 



47rAA 



122:12 I 



1 



-{^12) 



7rxi2 [1 - erf(a;i2)]| , 



(18) 



where X12 = |ri2|/AAi2. Notice that the color Kelbg potential approaches the color Coulomb potential at distances 
larger than AA12. What is of prime importance, the color Kelbg potential is finite at zero distance, thus removing 
in a natural way the classical divergences and making any artificial cut-offs, often applied (see, e.g., Ref. [isj). 
obsolete. This color potential is a straightforward generalization of the corresponding potential of electromagnetic 
Coulomb plasmas [2ll |. The off-diagonal color Kelbg potential can be approximated by the diagonal ones by means 
of «>i2(ri2, r'i2, Oi, Q2) « [*i2(ri2, ri2, Qi, Q2) + $i2(r'i2, r'12, Qi, Q2)]/2. 

Unfortunately such rigorous consideration of multiparticle density matrix for particles interacting by potentials 
unbounded below is not available. Therefore, following the experience gained in electromagnetic Coulomb plasmas, 
we use the following widely used ansatz [4^ \^ , which generalizes Eq. (|17p : 



p(') = p« exp 



1 ^ 

a4 e 



[■ ■ r ■ 



, Qi, Qj 



(19) 



Now we are able to construct R of Eq. ([T^ . The density matrix of noninteracting particles should be expressed 
in terms of determinants and permanents of single-particle density matrices in the standard way. Generalizing the 
electrodynamic plasma results (23j to the quark-gluon plasma case, we write approximate R 



R 



(r(0),rW,...r("+i);Q;{iV};/3 



N 



,r-',A/3 



exp{-/3c/}E nn- 

CT L/=i i=l 

per ||(/)(r("),r(o),A/3)||^ det ||0 (r'"' , rt") , A/3) | 



No 



a: 



3{n+l)Ng 



(A/3) 



det ||0(r("),r(o),A/3)| 



N„ 



A^("+'^^'(A/3) 



(20) 



In Eq. ((20|) the effective total color interaction energy is 



U = 



-1 n+l N 

+1^2 ^ , 



.(0 



, Qi, Qj 



'=1 hj{i¥=j) 

Other quantities in Eq. ([20]) are defined as follows: 



(21) 



(22) 



with Aa = \/27r/3/ma being a thermal wavelength of an a type quasiparticle (a = q,q,g). The antisymmetrization 
and symmetrization are taken into account by the symbols "det" and "per" denoting the determinant and permanent. 
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respectively. Eq. (|20|) is exact in the limit of n — > oo. Indeed, since each factor in Eq. p2|) has an error of order 
of l/(n + 1)^ with 6 > 1, the error of the whole product in the limit of n — > oo equals zero. Matrix </> (r, r', A/3) is 
defined by its matrix elements 

(r, r', A/3) = r^lMllI^^^ [6,^,^ + + Sa„,) Sf^j^ S.,^.^] (23) 
[zij [r, r', A^)] 

with 

(r, r' , A/3) = A/3 m, y^l + | r, - | '/ A/32 (24) 

defined in terms of the modified Bessel function A'2. These matrix elements are nonzero only for particles of the same 
type, i.e. ~ aj. Additional Kronecker symbols in spin, ai, and flavor, fi, indices of the particles are applicable only 
to quark and antiquark matrix elements. They prevent Pauli blocking for particles with different spins and flavors. 
The quantity describes the relativistic measure of trajectories in the color path integral. This measure is associated 
with relativistic operator of kinetic energy in Eq. In the limit of large mass this measure coincides with the 

Gaussian one used in Feynman- Wiener path integrals. Due to factors Sai.a ^fij- the matrix has a block structure 
corresponding to different types of particles and different flavors of quarks and antiquarks. Subscripts Na near det and 
per operations precisely refer to the corresponding blocks, which in case of quarks and antiquarks arc still subdivided 
into sub-blocks related to flavors. 

The dominant contribution to the partition function comes from configurations in which the "size " of the quasipar- 
ticle cloud of 'beads' is of the order of the Compton wavelength Ac — ^/rrii. Thus, this path integral representation 
takes into account quantum uncertainty of the quasiparticle position. In the limit of a large mass the spatial quasi- 
particle extension becomes much smaller than the average interparticle distance. This makes possible an analytical 
integration over the 'beads' positions by the method of steepest decent. As result the partition function is reduced to 
its classical limit involving point-like quasiparticles. 



III. WIGNER DYNAMICS 



We are going to use Wigner formulation of quantum mechanics for consideration of QGP kinetic properties. Let 
us review the underlying ideas of the Wigner dynamics for the simplest case, i.e. for nonrelativistic colorless system 
of particles [l^. The basis of our consideration is the Wigner representation of the von Neumann equation - a 
Wigner-Liouville equation (WLE). To derive the WLE for the simplest density matrix p{r,r',t) = 5'(r, i)^*(r', t) of 
the 3D N-particle system with ^ being an cigenfunction of a Hamiltonian operator H = "^i^i Pi /'m' + U, we introduce 
center-of-mass and relative coordinates in a standard manner: q = (r + r')/2 and ^ = r' — r. Note that all these 
quantities are 3iV-dimensional vectors. A Wigner distribution function (WF) is defined as 



w{p,q,t) 



1 



(2^) 



6N 



(25) 



Here and below products of vector quantities like are understood as scalar products of 3/V dimensional vectors. 
Using this definition, one can derive the WLE for w{p,q,t) [isl - isoj . Applying time derivative to definition and 
taking into account that 



z-f (r, t) = H^{r, t), i^**('^, t) = -H^*{r, t) 



we arrive at 

dw {p, q, t) 
dt 



1 



(27r)6JV 
1 



d^eMipO\ [H[T,t)-H[r',t) p{r,r',t) 

( i 

d^exp(ip^) <^ -— + - 



1 



(27r)6JV 

By means of integration by parts the first term in the braces can be transformed as follows 
1 



(2^) 



6N 



''^'^^^'P^^ ( j + 2' - 2' V = -m^^ 



(26) 



p[q+\,q-\,t 



(27) 



(28) 
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while for the second one we obtain the following expression 



1 



(27r)6W 
4 

6N 



exp(ip^) 



U q 



U q + 



dsw {p — s,q,t) / dq' U {q — q') sm{2sq') 



{2n) 

which results from substitution of p expressed in the form of inverse transformation to formula (^5 
This way we arrive at the following form of the WLE 



where 



dw p dw dU{q) dw 
dt m dq dq dp 

dU{q) dS{s) 



w(s,q) 



dq ds 



(27r) 



6N 



dsw{p — s, q, t) U) (s, q) , 



dq' U {q- q')sin{2sq') 



(29) 



(30) 



(31) 



In the classical limit, ^ — > 0, the r.h.s. of Eq. (|30p disappears and Eq. (PH]) is reduced to the classical Liouville 
equation. This is the reason why we extracted the term dU{q)/dq from the r.h.s. of Eq. ([50]) . 



A. Wigner Dynamics of Color Particles 



Let us consider dynamics of QGP quasiparticles additionally characterized by color variables Q and derive the 
color WLE for a density matrix p{r,r' ,Q,t) of the 3D N-particle system, where, as before, Q denotes color degrees 
of freedom of all quarks, antiquarks and gluons. Since color charges Q are treated classically, we consider only 
diagonal density matrix with respect to colors. Lideed, the Q variable already includes both canonical coordinate and 
momentum corresponding to classical color dynamics (see Appendix I). Therefore, the density matrix takes the form 

p(r, r', g, t) = p(r, r', Q, t) J] ^ iQ^ ~ ft W) (32) 

i 

where the product runs over all particles in the system and Qi{t) is a solution of the classical equation of motion 
for color Here p{r,r' ,Q,t) = \['(r, Q, t)\E'*(r', Q, t) is a quantum part of the density matrix with \1/ being an 
eigcnfunction of the Hamiltonian operator described by Eq. ([l} and Q arc already fixed c-numbers. 
Now the definition of the corresponding WF reads 

nj{p,q,Q,t) = j^^ J p(^q-^,q+i,Q,?je^Pid^, (33) 

The quasiparticles are also characterized by spin and flavor, which we do not explicitly include in the list of quasiparticle 
degrees of freedom. Notice that color degrees of freedom are also in the Wigner representation, since Q includes both 
color canonical coordinated and momenta, see Appendix L Now the WLE is defined by equation of the form [13, [Ml : 

dw dw ^dw ^^dw f , , ^ , 

^+v-+F- + T-^ J dswip-s,q,Q,t).(s.q), (34) 

where v = {v^} is 3N-dimensional vector of velocities of all quasiparticles, cf. Eq. ([3]), F = —dU^{q,Q)/dq is a set 
of the color-electric forces experienced by all quasiparticles, cf. Eq. T — {Tf} is an 8N-dimcnsional vector of 
driving forces in Wong's equation of motion for the color charge (jH), and uj (s, q) is defined by Eq. pip . 

The classical part of WLE (|34p . i.e. the l.h.s. of it, can be easily derived e.g. from the Wong's equations of 
motion (HJ-Q for the color-charged particles (see Ref. [51|). In particular, the term Tdw/dQ natirally results from 
{dQ{t) / dt)dw / dQ and Wong's equation Since we confine ourselves to classical dynamics of color, we do not need 
any further (quantum) consideration for it. The quantum space dynamics [i.e. the r.h.s. of Eq. p4p ] is derived 
completely in the same way as it was described above, see Eqs. (P7|) -([5n |) . with minor complications due to relativistic 
kinematics. 
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B. Wigner representation of time correlation functions 



In computations of transport properties, like viscosity, our starting point is the general Kubo expression for the 
canonical ensemble-averaged operator [s^ 



CBA{t) = Z-^Tr \e-P"l3 e'"' A e-'"'^ , 



(35) 



where B and A are quantum operators of dynamic quantities under consideration and Z {N, V, T) = Tr {e-^^} IS 
the canonical partition function. Frequently a symmetric time-correlation function is also used jssj : 



1. 



(36) 



where = t — i/3/2 is a complex- valued quantity including the inverse temperature j3 — 1/T. The Fourier transforms 
of CBA{t) and CBA{t) are related as [H^ 



Cba{^) = exp ( ) C'ba(w). 



(37) 



As a consequence, transport coefficients described by zero-frequency (w = 0) Fourier components can be obtained 
from the symmetric time-correlation functions, which may offer certain computational advantages. This symmetric 
form is used below. 

The Wigner representation of the time-correlation function in a 6-/V-dimensional space can be written as 



CBA{t) = (2^)-6W / dpqtiQdi^B{pqQ)A{pqQ)w(pqQ;pqQ;t;p'^ 



(38) 



where we introduced a short-hand notation for phase space points in (6N-|-8N)-dimesional space, pqQ and pqQ, with 
p, q and Q comprising the momenta, coordinates and color variables, respectively, of all particles of the system. Here 
A{pqQ) denotes the Weyl's symbol [4^ of the operator A: 



MpqQ) = / c^C exp(-iK) - |, Q 



A 



(39) 



and similarly for the operator S, while W [pqQ;pqQ:t; /3) is the spectral density expressed as 



W 



(pqQ ] pqQ -,1] 13 



= Z 



iHt* 



q--^.Q 



iHtc 



(40) 



In Eq. (|38p we silently assumed that operators A and B do not depend on spin variables. Therefore, summation over 
spins a and a can be safely moved to the definition of W. Here and below we do not explicitly write spin variables, if 
they are not essential. The time-correlation function CBA{t) is a linear functional of the spectral density W. Thus, 
the problem of its treatment is reduced to the consideration of the spectral density evolution. 

As it follows from Eqs. (|34p and Ref. [H^], the following system of the WL integro-differential equations describe 
the time evolution of the color spectral density W: 



dW _dW dW dW _ 

— +v—+F— + T-^ = I dsW [p- s,qQ]pqQ]t]P]uj{s,q,) , 

' + " I dsW{pqQ;p-s,qQ;t;l3]w{s,q), 



dt 



(41) 
(42) 



where as before u){s,q) is defined by Eq. pi[) . This equations are derived precisely in the same way as those of 
Eqs. (|34|) and ([27]) - ([30]) . only the Hamiltonian H appears here as a result of time derivation of exponent functions, 
e'^^ and e~*^*, rather than from application of equations of motion ([25)1 in Eq. (P7|) . Notice that while Eq. (PT|) 
describes evolution in the positive time direction, Eq. ([1^ specifies propagation in the reverse time direction. This 
happens because of the presence the direct time e~'^* and reverse time e*^* evolution operators in definition of the 
time-correlation function p6p . 
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Now using Eqs. (|^T|) and we can obtain an integral equation jisl - lsol . [s^] for W {pqQ\pqQ', t; {3 



dpoqofJ-QodpoqafiQoG {pqQ,pqQ,t;poqoQo,PoqoQo,Oj W{pQqoQo;poqoQo;t = 0,/3) 
W{p' - s,WQ';7q^': t'-: (3) u;{s,q') ~ WiWTW^P - s, qV'; t'; P) (43) 



with Green function 



G[pqQ,pqQ,t;p'q'Q',p'q'Q',t"j 
= 5{p- p{t; Wo', t')) S {q - q{t; WW, t')) 6 (Q - Q{t; WW, t')) 

X 5 (p-p(t;pW,iO) '^(9-9(^;pW^O)'5(Q-Q(^;pW (44) 

describing propagation of the spectral density along classical trajectories in positive time direction 



dp{t-p'q'Q\t') 1 
dt = 2^^'^^*^ 



dq{t\p'q'Q'.t') _ 1 



dt 2 



-v[p{t;p'q'Q',t% 



and in the reverse time direction 



dp{t;p'q'Q',t') Ip.^x 



dt 2 

v[pit;p'q'Q',t% 



dq{t-p'q'Q',t') _ 1^,^,,. , 



^ dt 2 
'MM^ = -iT(,0,), 

where {qQt) = \q{t]p'q'Q' ,t'),Q{t;p'q'Q' ,t')] and similarly for bared quantities. These equations of motion are 
supplemented by initial conditions at time t = Q 



p{t;poqoQo,0) = Po, ?(*; Po9oQo, 0) = go, Qit;PoqoQo,0) = Qo, (47) 

p{t;poqoQo,0) = Po, q{t'^PoqoQo,0) = qo, Q{t;poqoQo,0) = Qo- (48) 

and by initial conditions at time t ^ t' 



p{t'-p'q'Q',t')=p', q{t'-p'q'Q',t')^q\ Q{t' -p' q'Q' ,t') = Q , (49) 

pit';iwQ',t')^p', g(t';pW,i') Q(i';^W, *') = Q'- (50) 

In fact, Eqs. (|45p are Wong's equations of motion but written for half-time (t/2). Similarly, Eqs. (|46)) are half-time 
Wong's equations of motion reversed in time. This happens because the time correlation is taken between instants in 
the past and the future with the initial conditions fixed in between these instants, i.e. at t = the spectral density 
W^(w*3o' P^Qo' ^ = 0, /3) = Wo{pqQQ;pqQQ; f3) is fixed as it is described in the next subsection. 

Solution of the integral equation (|43p can be obtained in a form of iterative series. In this work, we take into 
account only the first term of this iteration series: 



w(pqQ;pqQ;t;(3) 

dpoqo/J'QodpoqoiiQoG (^pqQ,pqQ,t;poqoQo,poqoQo,0j WQ{pqQo;pqQo] P). (51) 
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Notice that if the initial Wo{pqQQ;pqQQ, P) is chosen appropriately [i^, i.e. such that it contains all powers of the 
Planck's constant, then the first term of the iterative series, i.e. Eq. (|5ip . describes propagation of a quantum initial 
spectral density along classical trajectories. Other (higher) terms describe propagation of the initial spectral density 
along the analogous trajectories but perturbated by momentum jumps resulted from the convolution structure of 
integral term in Eq. (143 p. From physical point of view these jumps relate to quantum ([p, q\) uncertainty (for details 
see discussion in Ref. j49l|). As it was found in Refs. [i^ [s^, UBIj the main contribution to WF comes from the 
trajectories without jumps, i.e. from the first term of the iterative series, if the motion takes place in a classically 
accessible region, which is the case here. Nevertheless calculations with higher-order iterative terms are in progress 
and will be reported elsewhere. 

Moreover, we replace the color Coulomb potential by the Kelbg one $ in quantities F and T defining propagation 
of the spectral density along classical trajectories, see Eqs. (P5|) and By this replacement we are able to take into 
account certain higher-order quantum terms of the iteration series presenting the solution of integral equation (j43p . 
From the practical point of view, it allows us to avoid problems due to singular character of the Coulomb potential. As 
known, the classical limit for multi-component Coulomb system does not exist since the stability of Coulomb systems 
is only provided by quantum effects. That is the physical reason for arising the color Kelbg potential in the patition 
function and Eq. (pT|) for initial condition (see below). 



C. Initial conditions 



The initial function Wq is expressed in terms of matrix elements of density matrix considered in section IIIBI 
Accordingly to Eq. (|40p for t = and definition of the density matrix p, cf. Eq. ([7]), we have 



X 5(Q-Q)fe,^. 



q-l,Q){q+l,Q 



(52) 



Thus, the problem is reduced to calculation of matrix elements of density matrix p, which is similar to that we did in 
sect. [TT] devoted to thermodynamics. Only now we need nondiagonal matrix elements rather than diagonal ones, as 
in sect. HIl 

As before, cf. Eq. (fT2|) . let us subdivide p(/?/2) into beards using the operator identity 



where the r.h.s. contains ti -I- 1 identical factors with A/3' = /3/[2(n -f 1)], so 



q+\.Q 



q-\,Q 



(53) 



where p^'^ (/ = 1, . . . , n -I- 1) are defined by Eqs. p3)) and ((T9| with 



= q 



and A/3 replaced by A/3'. The bar sign above p*-'^ means that these quantities depend on "bared" variables r^'^ 
Similarly, after antisymmetrization"' we obtain 



dF(i)...rfF("+i)p(i)p(2)p(3) ...^") 



p p~ p 



(54) 



* Symmetrization-antisymmetrization is performed only in one matrix element in Eq. H52|l . since it is enough to make the full expression 
under the trace II38I I to be symmetrized-antisymmetrized. 
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with 



^2 2 

where the "tilde" functions p^''^ depend on "tilde" variables. 

In Sect. [IT] we used the approximate expression for elements of the density matrix defined by Eq. p^ . Here it 
is not practical. An alternative approximation can be derived by means of the following expression for the density 
matrix operator 



exp « exp (^^^K^ exp (^-A^'U^^ exp jterms with [K,U'^]^ 



(55) 



compare with Eq. ([H)) . Notice that in fact exp{— (A/3'/2)iC} is a density matrix of the non-interacting system. This 
symmetrized form results in the following approximation to the matrix elements 



pM'-i),r('),Q;{iV}; A/372 



dqW p, (rC-i), Q; {N}; A/372) po {q^'\r^'\Q; {N}; Af3' 



X exp 



N 



(56) 



where is the diagonal part of the color Kclbg potential, see Eq. ^TE\\ . This approximate form has the same 
accuracy as that in Eq. (jl3l ) . At the same time it allows us to explicitly perform Fourier transforms in Eq. (j52p , 
since the ^ dependence now occurs only in the po factors. 

Thus, based on the above approximation for p*-') we are able to explicitly evaluate integrals over ^ and ^ 



-(1) 



(57) 
(58) 



where on the r.h.s. of these equations and below the marginal coordinates take already the following values 



and the complex- valued function ip is defined as 

N 



(p;r',r") = n(2^n exp 



r ■ — r 



(59) 



with A'; = yJirAfi'/rrii being the i-particle thermal wave lengths related to temperature 2/ A/3'. 
Substituting these expressions into Eq. ([5^ . we arrive at [46. ,54. ,55] 

Wo(mQ; pqQ\ 1^) = ^ J ^^^^ • • ■ '^^"^ ^^^^ • • ■ ^Q^"^^'' d^^^ ■ ■ ■ d^^"^ dq^^^ ■ ■ ■ dq^"+^^ 



^ {pqQ]pqQ-M^\..., 



tn).;i<i) ;T(n+i). 



(60) 



with 



* {pqQ;pqQ-M'\ . . . ,r("); g«, . . ., g^^+D . . .,7^^^;q^'\ . . . ,g<"+i);/3) 
[-13(U + U)]6(Q-Q) 



exp 
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'n+l 



E n^o(^^'"'^9^'^'Q;W;^/5'/2) Po(9^'\r«,Q;{iV}; A/372 

n 

n Po {^'-'^ > {^}; A/372) Po , , Q; {N}; A/372' 

E E E(-l)''"'^""^^?^^^s^o Q; {A^}; A/372) 5^,^ 5(5, P,P^P,5') 1^,.^ 



Pq Pq Pg 



where 



n+l TV 



2n 

1 1 

2n4 



n+l W 



lE E <i>.(qf'-qf,qf)-qy\Q. 



Applying the notation of Sect. II (cf. Eqs. and Eq. (I^^l) ). wc finally arrive at 

* (^;p^;r(i), . . . ,r("); g^^^, . . . ,g("+i';?<i), . . . • • • ,g^"+'^;/3) 

= exp{-/3(c7+J7)} (5(0 -Q) 



n+l / N 

n H' 

cr,cr L^ — 1 \i— 1 

TV 



E 



r('-i) 77' 



g«, A/372) ^ (^n^,, (g«,r«, A/372) ^ 



n n'^»K^'"'^5<'\ A/372) n'^^'(9^'^^'^ A/372) 

,Z=1 \i=l / \i=l 

det ||(/.(g<"+i),?<"+i),A^72)||^ det ||(/) A^72) ||^ 



A 



6(n+l)7V, 



(A/372) 



A^("+^''^7A/372) 



per ||(/) (5^' 



+ 1) p.ri+1) AO/ 



A/372) I 



A' 



6(n+l)iVg 



(A/372) 



with 



(61) 



(62) 
(63) 



(64) 



In the limit n — ?> 00 this expression exactly gives the product of matrix elements in Eq. ([5^ in the form of path 
integrals multiplied by a limiting expression of the ip functions. According to the Lebesque-Dirac delta theorem (see 
Appendix II), the p'-^p^^'ip products in integral in the limit n — 00 are equivalent to the to the following real 
valued expressions 



{r. 
(p. 



TV 



n-4-oo I 



i=l 

N 



n((2A?)^/^xp 



i=l 



n((2An'/'cxp 



r fi2\2i 
Pi \ 


)(^ 


27r _ 




Pi \ 




27r _ 





(65) 
(66) 



where A,; = 7r/3/2mi. Analytical integration over the delta function simplifies the final path integral used for further 
computation of the real-valued Wq by means of the Monte-Carlo method. 
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IV. SIMULATIONS OF QGP 



The developed approach is apphed to the QGP at zero baryon density = 0). Then assumption on equal quark 
masses (see point IV in subsect. Ill A[) immediately implies equal fractions of quarks and antiqiarks of different flavors: 
iVtj = Nfi = Ng = Nq/3 = Nu = Nd = Ns = Nq/3. Ideally the parameters of the model should be deduced from 
the QCD lattice data. However, presently this task is still quite ambiguous. Therefore^Jn the present simulations we 
take only a possible set of parameters. Wc use so called "one-loop analytic coupling" [5^, [131 



11 - i2/3)Nf 



1 A2 



QCD 



HQV^Icd) ^%cd - 



(67) 



where Q is the momentum transfer, Aq^d = 206 MeV is the QCD scale and N f = 3 is the number of flavors. The 
analytically generated non-perturbative contribution ^^cdI{^%cd ~ Q^) subtracts the unphysical Landau pole in 
a minimal way, yielding a ghost-free behavior which avoids any adjustable parameter. This coupling agrees with a 
great body of experimental data [H, [s^- As it is usually done in thermal field models, we substitute Q by 2ttT 
to use this coupling in our simulations. The resulting as{T) is displayed in the left panel of Fig. [1] and compared 
with QCD-lattice coupling deduced from a short-distance behavior of the singlet free energy (ssj and from spectral 
density of heavy-quark correlator [59| . As seen, the running coupling deduced from experimental data is close to 
those obtained in the lattice QCD. Notice that determination of as{T) in lattice QCD simulations is quite indirect. 
Therefore, different indirect methods naturally give somewhat different results. 




FIG. 1: Input data for the PIMC simulations. 

Left panel: Running coupling constant versus temperature fitted to experimental data [bgI. [STI (solid line). Points present the 
coupling deduced from lattice QCD simulations in Refs. [H^ [lattice (Bielefeld 2012)] and [Hy] [lattice (Banerjee et al. 2012)]. 
Right panel: Mass-to-temperature ratio for quark and gluon quasiparticles versus temperature. Points are values used in 
simulations. The solid lines are smooth interpolations between points. 

The quasiparticle masses were chosen to reproduce the pressure obtained in lattice QCD calculations 0, Hj. The 
T-dependence of these masses is pres ented in Fig. [Ij right panel). When choosing masses we kep t in mind constraints 
resulting from lattice QCD data |35l . IstI . |60| and QCD- motivated quasiparticle models [6l|, HI]. While gluon masses 
used in this paper well comply with those deduced from both lattice QCD data [HjliOl and quasiparticle models, this 
is not the case for quark masses. Our quark masses agrees with values required for quasiparticle fits j6ll l62j of the 
lattice thermodynamic properties of the QGP: iriq/T ~ 1.5 2.5. At the same time they are appreciably lower than 
those in old lattice data [s^: rUq/T ~ 4, and higher than iriq/T ~ 0.8 reported in newer lattice calculations (37j . 



A. Equilibrium Properties 



Figured] (left panel) demonstrates the quality of reproduction of the equation of state (EOS), i.e. the pressure 
versus temperature, achieved with the above discussed input data. The reference EoS (filled points in the left panel 
of Fig. is taken from QCD lattice simulations of the QGP The quality of the reproduction obviously depends 
on the degree of accurate tuning of the input data, i.e. the quasiparticle masses. However, the PIMC scheme itself 
also produces certain errors. If there are metastable states of the system, convergence of calculations becomes poor 
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because of jumps between stable and metastable states. This is a typical situation when the system approaches to a 
point (or a range) of a phase transition. Precisely this happens at the lower end of considered temperature range. The 
shaded aria in the left panel of Fig. [2] indicates these uncertainties of the PIMC calculations. Figure [2] also presents 
the entropy S/T^ and trace anomaly (e — 3P)/T'^ of the QGP. These quantities are calculated accordingly to Eqs. 
((8l)- pT|) . In order to avoid the numeric noise, the derivative of a smooth interpolation between the PIMC points 
(solid line in the left panel of Fig. [2]) was taken. Though agreement with the lattice data looks worse for the entropy 
and trace anomaly, in fact it is the same as that for pressure. Differentiation operations in Eqs. ^ and (|lip make 
differences between PIMC results and lattice data more pronounced. 




100 200 300 



400 500 
T[MeV] 



100 200 300 400 500 
T [MeV] 



100 200 300 400 500 
T [MeV] 



FIG. 2: Pressure (left panel), entropy (midle panel) and trace anomaly (right panel) scaled with corresponding powers of 
temperature versus temperature from PIMC simulations (open squares). These are compared with lattice data of Refs. 0, [H| 
(filled circles). The solid line in the left panel is a smooth interpolation between the PIMC points. Entropy and trace anomaly 
are calculated based on this smooth interpolation. The shaded aria in the left panel indicates uncertainties of the PIMC 
calculations. 



Having calibrated the model by reproducing the EoS, we can proceed to predictions. First, let us consider internal 
properties of the system. To characterize physical conditions and interplay of interaction and degeneracy in Fig. [3] a 
degeneracy parameter Xu for 'up' quarks and a plasma coupling parameter F are presented: 

Xu = n^\l r = £^, (68) 

where the thermal wave length A„ was defined in the previous sectios (see text after Eq. (|22[) ). Uu is density of u 
quarks, = 3/(47rn) is Wigner-Seitz radius, n is the density of all quasiparticles (quarks, antiquarks and gluons), 
and ^2 is the quadratic Casimir value averaged over quarks, antiquarks and gluons, ^2 = — 1 is a good estimate 
for this quantity. The plasma coupling parameter is a measure of ratio of the average potential to the average kinetic 
energy, and the degeneracy parameter Xu indicates wheather a system is classical (x„ <C I) or quantum (xu !)• It 
turns out that F and Xu are of order unity which indicates that the QGP is a strongly coupled quantum (xu > 
liquid (r ^ 1) rather than a gas. 

To clarify interplay of interaction and degeneracy let us consider spatial arrangement of the quasiparticles in the 
QGP by studying a pair distribution function (PDF) gab{R) defined as 

g,b(|Ri-R2|) = C^] E Sa„aSa„b^JdrdtiQSiRi~r,)5{R2-r,)pir,Q,a;{N};P), (69) 

where at and Uj are types of the particles (= q,q or g). The PDF gives a probability density to find a pair of particles 
of types a and 6 at a certain distance R = |Ri — R2I from each other. The PDF depends only on the difference of 
coordinates because of the translational invariance of the system. In a non- interacting classical system, gab{R) ^ I, 
whereas interactions and quantum statistics result in a redistribution of the particles. At temperatures T = 525 MeV 
and T = 193 MeV the PDF averaged over the quasiparticle spin, colors and flavors are shown in Fig. 21 

At distances R > 0.2 or 0.3 fm, depending on the temperature, all PDF's are practically equal to unity (Fig. U]) like 
in ideal gas due to the screening of the color Coulomb interaction. A drastic difference between qq and gg PDF's (the 
qq PDF is identical to the qq one) occurs at short distances. Here the gluon-gluon and gluon-quark PDF's increase 
monotonically when the distance goes to zero while the qq and qq ones remain uncorrclatcd. One of the physical 
reasons of the PDF difference is spatial quantum uncertainty and different properties of Bosc and Fermi statistics of 
gluon and (anti)quark quasiparticles. Uncertainity in particle localization is defined by the ratio T/m. Localization 
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FIG. 3: Quark degeneracy parameter Xu and the plasma coupling parameter F [see Eq. (|68p] versus temperature. Lines are 
smooth interpolations between points. 
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FIG. 4: Pair correlation functions of identical (left panels) and different (right panels) quasiparticles at temperatures T 
525 MeV (upper panels) and T — 193 MeV (lower panels). 



is better for heavier gluon quasipartcles. Fermi statistics results in effective quark-quark and aniquark-antiquark 
repulsion, while Bose one results in effective qluon-qluon attraction. Oscillations of the PDF at very small distances 
R < 0.1 fm are related to Monte-Carlo statistical error, as probability of quasiparticles being at short distances quickly 
decreases. 

However, the qq and qq pair correlation functions reveal practical absence of fermi repulsion. This happens because 
another physical reason comes into play. Strong interaction between quasiparticles reduces the influence of the 
degeneracy in the region of Xu ~ 1- This interaction is dominated by attraction at short distances. Indeed, the QGP 
lowers its total energy by minimizing the color Coulomb interaction energy via a spontaneous "anti- ferromagnetic"- 
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FIG. 5: Gluon-gluon and quark-antiquark pair correlation functions multiplied by i? at T = 193 MeV. 



like ordering of color vectors, i.e. the color vectors of nearest neighbor quasiparticles become anti-parallel. Similar 
absence of fermi repulsion was observed in hydrogen plasma at x ^ ^ [47| . This short-distance attraction is stronger for 
gluon-gluon and gluon-(anti)quark pairs than for (anti)quark-(anti)quark ones because of the corresponding difference 
in values of quadratic Casimir invariants q2 (see Appendix I), which determine the maximal values of the effective 
color charge products \Qi ■ Qj\ in color Kelbg (Coulomb) potentials: For gluon-gluon pairs \Qg ■ Qgl^^^^ = 24, for 
gluon-(anti)quark pairs \Qg-Qq\^ax = IQg'Qqlmax ~ ^'^'^ ^^"^ (anti) quark- (anti) quark pairs \Qg-Qg\^^^ = 
\Qq ■ Qqljnax ~ 1^9 ' Q'^\max ~ Stronger gg attraction additionally enhances correlation of the gluon-gluon pairs 
at short distances. At the same time the short-distance attraction is the only reason of the gluon-(anti)quark short- 
distance correlation. 

The short-distance correlation implies formation of the gluon-gluon and giuon-(anti)quark clusters, which are uni- 
formly distributed in space. In case of the gluon-gluon clusters we can even speak about gg bound states, i.e. glueballs, 
as it is seen from Fig. [S] The product R^gabiR) is proportional (up to constant factor) to a probability to find a 
pair of quasiparticles at a distance R from each other. As known from consideration of hydrogen plasma [47|, a max- 
imum in ggg{R) signals population of a bound state. For comparison, the quark-antiquark correlation function, 
i.e. R^ gqg{R), is also presented in Fig. [S] It demonstrates that there are no bound meson-like states. We can only 
speak about weak meson-like clustering at lower temperatures, see short-distance qq correlation at T=193 MeV in 
Fig. m Possible existence of medium- modified meson-like bound states was actively discussed some time ago, e.g ., 
in Ref. [gJ] and later in Refs. [H, llBl based on results from lattice QCD calculations of spectral functions j67l [68j. 
Our result supports conclusion of Ref. [6^ on the absence of qq bound states above the temperature of the phase 
transition. This finding is in contrast to our previous results on SU(2) group [30l - l33| . There well pronounced bound 
qq states were found just above the critical temperature, which however quickly dissolved with the temperature rise. 
This happens because the SU(3) plasma turns out to be essentially denser than the SU(2) one, which is a consequence 
of a stronger effective attraction between constituents. As a result, possible bound states in the SU(3) plasma just 
melt. To verify the relevance of all above discussed trends, a more refined color-, flavor-, spin-resolving analysis of 
the PDF's is necessary. This work is presently in progress. 

Details of our path integral Monte-Carlo simulations have been discussed elsewhere in a number of papers and 
review articles, see, e.g. Refs. j47l . [63| and references therein. The main idea of the simulations consists in con- 
structing a Markov chain of configurations which differ by the quasiparticle coordinates. Additionally to the case of 
electrodynamic plasmas, here we randomly sample the color quasiparticle variables according to the SU(3) group Haar 
measure with the quadratic and cubic Casimir conditions until a full convergence is achieved. Another important 
difference from the case of the electrodynamic plasmas consists in using the relativistic measure in Feynman- Wiener 
path integrals associated with relativistic kinetic energy operator instead of the conventional Gaussian one arising 
from the non-relativistic operator. For simulations we use a cubic simulation box with periodic boundary conditions. 

At the first step a dominant, i.e. maximal, {A^}-term in the sum of Eq. ^ was determined by calculations in 
grand canonical ensemble. This term is indeed dominant in the thermodynamic limit of the box volume — > oo. Then 
all the calculation was performed for a canonical ensemble with thus determined {A^}. In practical calculations we 



always take the number of particles TV = Nq 



Ng ~ 126 and adjust the above determined proper densities 



of species by varying the total volume V of the box. The number of beads n = 20 for each particle was used. On 
one hand, errors in MC calculations of thermodynam ic q uantities related to the finite number of particle in a system 
with periodic boundary conditions are of order 1/A^ |47| . On other hand too large number of particle presented by 
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large number of beads requires large computer resorces. Our choice of particle and beads numbers is a compromise 
between acceptable accuracy and available computer resorces. Test variation of the number of beads from 15 up to 
50 practically does not change results. 



B. Transport coefficients 



An important aspect of the strongly coupled QGP is its transport properties which strongly differ from those we 
would expect for weakly coupled plasmas. We use the developed approach based on Wigner formulation of quantum 
mechanics to calculate the QGP transport properties at strong coupling. In particular, we calculate the QGP self- 
diffusion constant and shear viscosity, as these quantities can be compared to respective values deduced from analysis 
of experimental data heavy-ion collisions and also predictions of the lattice QCD computations. More precisely, 
summary of shear viscosity deduced from analysis of experimental elliptic flow is presented in Ref. (73| , in Ref. (74 | 
an extensive review of theoretical works on viscousity is done, while the heavy-quark diffusion constant is available 
from the experiment analysis [7l|, [T^I and QCD lattice computations [H^, [t^]- We anticipate that the self-diffusion 
and heavy-quark diffusion constants are compareable by the order of magnitude. 

A natural way to obtain these transport coefficients is use of the quantum Green-Kubo relations. These relations 
give the transport coefficients in terms of integrals of equilibrium time-depended correlation functions. According to 
the Eq. psp a self-diffusion constant D is the integral of the velocity autocorrelation function 

D lim D{t), 
m = ^ j dT{v{0) ■ vir)) 
(i>(0)-«(r)) = {2^)-'''' j d^di^wi^-^-T-fi) v{p)-v{p), (70) 
where the product of 3-velocities is 

N _ ~ 

v{p)-v{p)^y. _ ^'' ^1 — (71) 

Figure [5] shows examples of the velocity- velocity autocorrelation and its antiderivative function. The self-diffusion 
constant is a limiting value of the related antiderivative function at t oo. 
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FIG. 6: Velocity autocorrelation function (left panel) and temperature-scaled self-diffusion function T D(t) (right panel) versus 
time for two temperatures. 



Analogously, the Green-Kubo relation for the shear viscosity is the integral of the autocorrelation function of the 
stress-energy tensor 

7y = lim ri{t), 

1 
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{(JxyiO) (Jxy{T)) = {2t:Y^^ J dpqfiQ dpqfiQ W [pqQ] pqQ] r; (3^ a^y (pqQ) (Txy (pqQ^ , (72) 

where the off-diagonal stress-energy tensor is 

, ^^ Pt,x Pt^y dU jqQ) 

a,y{pqQ)^^^===--^q,,^x^-—, (73) 

here qij = Qi — q^, and U is the sum of the color Kelbg potentials defined by Eqs. (|62|) and (|63| with n = 0. 
Examples of the stress-energy-tensor autocorrelation and its antiderivative functions arc presented in Figure [71 The 
shear viscosity is defined by limiting value of the related antiderivative function at t — > cx3. 





FIG. 8: Left panel: The self-diffusion constant [self-diffusion, Wigner dynamics] as a function of temperature, compared 
with heavy-quark diffusion constants predicted by QCD lattice computations [c-quark, lattice (Ding et al.2011)] and [H^ 
[c-quark, lattice (Banerjee et al. 2012)], and deduced from analysis of experimental data [tJ] [c-quark, "exp." (Gossiaux at al. 
2012)] and ^ [c-quark, "exp." (He at al. 2012)]. 

Right panel: The ratio of shear viscosity to entropy density as a function of temperature [shaded aria marked "Wigner 
dynamics"]. Horizontal dashed lines indicate the range of constraint on the viscosity-to-entropy ratio deduced from numerous 
hydrodynamical simulations of heavy-ion experimental data, as summarized in Ref. [t^ ]. 

The self-diffusion constant and the viscosity-to-entropy ratio are presented in Fig. [5] as a function of temperature. 
The entropy density s is taken from results of our PIMC calculations presented in Fig. [51 Our results are presented in 
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the form of bands. The width of the bands represents the theoretical uncertainty associated with the oscillations of the 
antiderivative functions at large times, see Figs. [n]and[71 Slowly decaying oscillations of the time correlation functions 
are inherent in liquid-like systems of strongly interacting particles in contrast to exponentially decaying oscillations 
in gas-like systems. As known from hydrogen plasma, these oscillations arise because of quasiclosed chaotic orbits 
and are caused by strong interparticle interaction. In liquids these oscillations decay according to a power law rather 
than exponentially. Therefore, extremely long (in time) trajectories are required for more accurate calculations of the 
diffusion constant and viscosity. Due to CPU time limitations, we had to stop our computations before the decay of 
these oscillations was completed. 

Unfortunately, the self-diffusion constant is unavailable from other calculations. Therefore, we compare it with 
heavy-quark diffusion constant, anticipating that these are of the same order of magnitude. The heavy-quark diffusion 
constant is available from recent QCD lattice computations [1^ [t^], which are presented in the left panel of Figs. HI 
Our results (labeled as "Wigner dynamics" ) well agree with lattice data of Ref . [s^ , while essentially overestimate 
those of Ref. [tO ] . The heavy-quark diffusion constant is also available from analysis of experiments on heavy-quark 
quenching in ultrarelativistic heavy ions collisions at RHIC. For comparison we took two recent works on the subject 
[tH [72l |. Here the results are also rather diverse. Estimates of Ref. [tH are well comparable with our result and 
lattice data of Ref. , while estimate of Ref. [7l|, [t^ is considerably lower and better conforms to lattice results of 
Ref. H. 

Our results on the shear viscosity are presented in the right panel of Figs. HI As seen, the obtained values of 
viscosity are in the range of those deduced from the analysis of experimental elliptic flow in ultrarelativistic heavy 
ions collisions at RHIC, as summarized in Ref. [t^. Lattice data on the shear viscosity in the realistic case of the 
SU(3) group are not available, an extensive review of theoretical works on viscosity within QCD-motivated models 
is done in Ref. (t^ . As seen, the minimum of the viscosity-to-entropy ratio is reached at a temperature above the 
expected phase transition rather than in the phase transition range, as is commonly expected. This minimum turns 
out to be quite shallow. The value of the viscosit y-to -entropy ratio at the minimum is very close (from above) to 
the lower bound of ij/S = l/in for this quantity [7a], often referred to as the KSS bound. With the temperature 
decrease, i.e. towards the hadronic phase, the viscosity rapidly rises. 



V. CONCLUSION 



In this paper we demonstrated that color quantum Monte-Carlo (PIMC) simulations based on the quasiparticle 
model of the QGP are able to reproduce the lattice equation of state at zero baryon chemical potential at realistic 
model parameters (i.e. quasiparticle masses and coupling constant) even near the critical temperature and also yields 
valuable insight into the internal structure of the QGP. In our simulations we have introduced a new relativistic path 
integral measure and have developed a procedure of sampling color quasiparticle variables according to the SU(3) group 
Haar measure with appropriate Casimir conditions. Unfortunately, convergence of our calculations becomes poor in 
the range of the expected phase transition because the scheme suffers from jumps between stable and mctastable 
states which turn out to be almost equally probable in this range. Our results indicate that the QGP reveals quantum 
liquid-like (rather than gas-like) properties up to the highest considered temperature of 525McV. 

Short-distance correlations in the computed pair distribution functions of gluon-gluon and gluon-(anti)quark pairs 
displays the formation of clusters. In case of the gluon-gluon clusters we can even speak about gluon-gluon bound 
states, i.e. glueballs, at temperatures just above the phase transition. The possible existence of medium-modified 
meson- like bound states was actively discussed some time ago [6^ [6^ . Our result supports the conclusion of Ref. 
[69] on the absence of qq bound states above the temperature of the phase transition. This finding is in contrast to 
our previous results on the SU(2) group [30l - [33j . There well pronounced bound qq states were found just above the 
critical temperature, which however quickly dissolved with the temperature rise. This happens because the SU(3) 
plasma turns out to be essentially denser than the SU(2) one, which is a consequence of a stronger effective attraction 
between the constituents. As a result, possible meson-like bound states in the SU(3) plasma just melt. 

The PIMC method is not able to yield transport properties of the QGP. A way to access these is to develop a classical 
color molecular dynamics simulation [l^ , where quantum effects are included phenomenologically via a short-range 
potential. In contrast to these classical MD simulations we have developed a more rigorous approach based on 
the combination of Feynman and Wigner formulations of quantum dynamics. The basic ideas of this approach have 
been briefly reported in Ref. [3^ . In this paper we gave a more detailed description. In particular, this approach 
allowed us to calculate the self-diffusion coefficient and the viscosity of the strongly coupled QGP. Since the self- 
diffusion constant is unavailable from other calculations, we compared it with the heavy-quark diffusion constant, 
anticipating that these are of the same order of magnitude. The heavy-quark diffusion constant is available from 
recent QCD lattice computations [sol . [70j and also from an analysis of the heavy-quark quenching in experiments on 
ultrarelativistic heavy ions collisions at RHIC. For comparison we took two recent works on such an analysis (tiIIt^. 
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Unfortunately the above mentioned lattice and heavy-quark-quenching results are rather diverse. Our self-diffusion 
constant well agrees with lattice data of Ref. [s^ and estimates of Ref . [72| , while essentially overestimates those of 
Refs. [zl and [zl. 

Our results on the shear viscosity are in the range of those deduced from the analysis of the experimental elliptic flow 
in ultrarelativistic heavy ions collisions at RHIC, as summarized in Ref. [73| . i.e. in terms of the viscosity-to-entropy 
ratio, 1/Att < rj/S < 2.5/47r, in the temperature range from 170 to 440 MeV. The minimum of the viscosity-to-entropy 
ratio is reached at a temperature 300 MeV) , above the expected phase transition rather than in the phase transition 
range as commonly expected. This minimum turns out to be very shallow. The value of the viscosity-to-entropy ratio 
at the minimum is very close (from above) to the lower bound of r]/S ~ for this quantity [75|, i.e. to the KSS 
bound. With the temperature decrease, i.e. towards the hadronic phase, the viscosity rapidly rises. 

Our present analysis is still confined only to the case of zero baryon chemical potential. Simulations at nonzero 
baryon chemical potentials are in progress. 

We acknowledge stimulating discussions with P. Levai, D. Blaschke, R. Bock, H. Stoecker, and D.N.Voskresensky. 
Y.B.I, was partially supported by grant of the Russian Ministry of Science and Education NS-215.2012.2. 



Appendix I: Integration over SU(3) group Haar measure 



In this appendix we explain details of integration over SU(3) Haar measure dfiQ in Eq. 
color charge in case of the SU(3) group is [ifl, l76l - [78j 

d^J.Q = (fQ 5{QaQa - q2) S{dabcQaQbQc - qs) 



|). The measure for single 



(74) 



with summation over a,b,c = 1, . . . , 8 and constants dabc given in Table 1. For the SU(N) group the quadratic q2 
and cubic 53 Casimirs are 92 = (-/V^ — 1)C^ with C2 = N for qluons and C2 = 1/2 for quarks and antiquarks, 93 = 
for gluons and 93 = {N"^ — 4)(7V^ — l)/4 for quarks, for antiquarks 93 has opposite sign. In fact, the normalization 
constant Cft depends on 92 and 93 Casimirs. 

For random sampling of the Q variable in Monte-Carlo integration in Eq. ^ we change to the related canonical 
Darboux variables for the SU(3) group. The set of the canonical variables [07r] [((/iq,, tTq., Q = 1,2,3] is defined by 
the canonical Poisson bracket 

_ dA dB dA dB dA dB dA dB 



{A,B}, 



yr 00 an dn d(j) 
where r and p are conventional coordinate and momentum, respectively, and obey 

The color variables Qa form a representation of SU(3). In terms canonical variables their Poison bracket reads 

where fabc are the structure constants of SU(3) given in Table 1. 

The explicit transformations to canonical variables are given [lol . [76l - [78| by expressions: 



(75) 

(76) 
(77) 



Qi = 7r^7r_ cos 01, (32 = 7r+7r_ sin(/)i, 

Qi = C++TT+A + C+-TT-B, Q5 = S++n+A + S+-7r-B, 



in which we have used definitions: 



C±± 

and A and B are given by 

A^ 

B = 



71'+ = V^s + Tn 
'1 



-TT,B, 



V302±, 



h = S- 



-A-S—TT+B, 



S±± — sin 



i (±01 + V3(j)2 ± 03) 



(78) 



(79) 




J1-J2 



7r3 



71-2 

x/3 



Ji + 2 J2 



7r3 



73 



2 Ji + J2 



T^3 



7^2 

V3 



J2- Jl 



T^3 



TT2_ 

V3 



Ji + 2 J2 



7r3 



7r2 
73 



2 Ji + J2 



7r3 



7r2 
73 



(80) 
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In this expression the set Qi,Q2,Q3 forms an SU(2) subgroup with quadratic Casimir Q\-\- Q\ + Q\ = TTg. Let us 
note that two Casimirs depend only on Ji and J2. They can be computed using the values given in Table 1 as 



J1J2 + ./, 



1 



dabcQaQbQc = (^1 - J2) (Jl + 2 J2) (2 Jl + J2) 
io 



(81) 



The phase space color measure for SU(3), given in Eq. (|74p can be transformed to the new coordinates through 
use of Eq. ([75)1 and evaluation of the Jacobian 



Then the measure reads 

djiQ ~ d(f>ic 



\/3 

i>^d-nidTi2dTT3dJidJ2 — JiJ2{Ji + J2) 
4o 



(82) 



1 



-^(Ji - J2)(Ji + 2J2)(2Ji + J2) - 93 
io 



(83) 



Since the two Casimirs are independent, the 5-functions fix both and J^. After integration over and J2 the 
Eq. ([55]) gives a proper canonical volume element dcjidi:. Thus applying Metropolis algorithm to (pi: variables we can 
construct Markovian chain in (pi: phase space and obtain random color variables Q for calculation partition function 
according to the SU(3) group Haar measure with two Casimir conditions. 



TABLE L Non-zero fabc and dabc constants of the group SU(3) 



fabc 


/123 


/l47 


/l56 


/246 


/257 


/345 


/367 


.Ass 


/678 


values 


1 


1 

2 


1 
2 


1 
2 


1 
2 


1 
2 


1 
2 


2 


V3 
2 



dabc 


diis 


dim 


dl57 


^228 


^247 


^256 


^338 


^344 


0(355 


1^366 


1^377 


d448 


^558 


^668 


^778 


^888 


values 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 
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1 


1 


1 




2 
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V3 


2 


2 


V3 


2 


2 


2 


2 


2^3 


2^3 


2^3 


2\/3 





Appendix II: Lebesque-Dirac delta theorem 

Let / be a summable function of real argument such that J f{x)dx = Iq. 
Then Mf{M{x - x')) IoS{x - x') when M 00. 
Proof: Let (p be any test function. Then 

lim I Mf{M{x-x'))ip{x)dx= lim [ f{s)(p{s/M + x')ds 

M— )-oo J AZ-i-oo J 

f{s){ lim ip{s/M + x'))ds^ip{x') [ f{s)ds^Io(p{x') 
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